% input the txt file's name to get the longitudinal, perpendicular and total correlation  
function [ ] = computeCorrelationExp( dirName,fileNum,dataMode )
global computeCorrC

[position,velocity,num]=readText(dirName, fileNum, dataMode); % use the function "readText"

S = sum(velocity) / num; % now S is the average velocity vector 
SL = norm(S);
normalS = S./norm(S);
distance = zeros(num,num);
   for i= 1 : num
        for j = 1: num
            distance(i,j) = norm(position(i,:)-position(j,:));
        end
    end % acquire distance and correlation between i and j
sL = velocity*normalS';
pai = velocity - sL*normalS;

[logicI,~] = distinguishBoundary(position,distance);
 %logicI = true(1,num);
numI = sum(logicI);
pai_I = pai(logicI,:);
sL_I = sL(logicI,:); 
distance_II = distance(logicI,logicI);% an important correction !

sL_approx = 1-0.5*(pai_I(:,1).^2+pai_I(:,2).^2+pai_I(:,3).^2);

disp('approximational error = ');
disp(sum(abs(sL_I-sL_approx))/length(sL_I));
disp('sL difference = ');
disp(sum(abs(sL_I-SL))/length(sL_I)); % two important error data, average? maximum?

corrP = pai_I * pai_I';
corrL = (sL_I-SL) * (sL_I-SL)';
% get the matrix of correlation for further computing
% ----------------------Compute Correlation------------------------%
maxR = floor(max(max(distance_II)))+1;
deltaR = maxR / 20;
maxN = floor(maxR/deltaR); 


CorL = zeros(1,maxN);
CorP = zeros(1,maxN);
count2 = zeros(1,maxN);
for i = 1 : numI
    for j = 1 : numI
        temp = floor((distance_II(i,j)/deltaR))+1;
        if temp < maxN
        CorP(temp) = CorP(temp) + corrP(i,j);
        CorL(temp) = CorL(temp) + corrL(i,j);
        count2(floor((distance_II(i,j)/deltaR))+1) = count2(floor((distance_II(i,j)/deltaR))+1)+1;
        end
    end
end
CorP = CorP./count2;
CorL = CorL./count2;


% Cor = CorL+CorP;
% --------------------------Plot-----------------------------%
rLabel = deltaR*(1:maxN);
numLabel = 1:10:numI;
if dataMode == 1
    figure(5)
    quiver3(position(:,1),position(:,2),position(:,3),velocity(:,1),velocity(:,2),velocity(:,3));
   
elseif dataMode == 2
    figure(5)
    quiver(position(:,1),position(:,2),velocity(:,1),velocity(:,2));
    axis equal
end
figure(2)
subplot(1,2,1)
plot(rLabel,CorL,'ro');
hold on
plot(rLabel,zeros(1,maxN));
title('longitudinal correlation of r')
subplot(1,2,2)
plot(rLabel,zeros(1,maxN));
hold on
plot(rLabel,CorP,'ro'); 
title('perpendicular correlation of r')

if computeCorrC ==1
    figure(3)
    subplot(1,2,1)
    plot(numLabel, pai_I(numLabel),'ro')
    hold on
end
% subplot(1,3,3)
% hold on
% plot(r,Cor,'ro');
% title('total correlation of r')  % draw the plot

% figure(4) % test
% hold on
% plot(sL,sL_approx,'r.');
% plot(sL,sL)
% hold off
    